source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/SLE_TRPM2_MfMo/'

kn.pal2 <- c('#d24040','#485ffc')

# zheng22 SLE skin -------
sobj <- read_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_derm_epi.rds')

sobj |>
  write_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_derm_epi.rds')

## big umap ----------
sobj |> DimPlot(group.by = 'manual_fine', cols = 'Paired', raster = F) +
  ggtitle('SLE Skin') +
  theme_jpub()

sobj |>
  ggplot(aes(umap_1, umap_2, fill = manual_fine)) +
  geom_bin2d(bins = 512) +
  theme_jpub(theme_classic) +
  scale_fill_brewer(palette = 'Paired') +
  labs(title = 'SLE skin', fill = '')
  
publish_pdf('SLE.skin.umap.pdf', width = 70)

sobj |> DotPlot(pbmc_markers, cols = 'RdYlBu', group.by = 'seurat_clusters')

sobj |> DimPlot(group.by = 'orig.ident')

sobj |> DimPlot(group.by = 'tissue')

sobj |>
  DotPlot(c('TRPM2','CD14'), group.by = 'manual_main',
          cols = 'RdBu') +
  labs(x = 'Gene', y = 'Cell type', title = 'SLE skin')

## FeaturePlot TRPM2 SLE/HC ---------
sobj |> FeaturePlot('TRPM2', cols = c('lightgrey', 'red'), order = T,
                    split.by = 'sle', raster = F)

g0 <- last_plot()

g0 +
  plot_annotation(title = 'TRPM2 expression') & theme_jpub & NoLegend()

publish_pdf('sle.skin.trpm2.featureplot.pdf', width = 100, height = 50)

## DC(mf) in epidermis -------
sobj <- read_rds('mission/FPP/xiangya_sle_scRNA/mix3.epi.sle.rds')

sobj.dc <- sobj |> filter(manual.main == 'DC')

sobj.dc <- sobj.dc |>
  quick_process_seurat(batch = c('orig.ident','batch'), skip_norm = T)

sobj.dc |>
  DimPlot(cols = paired.pal, label = T, label.box = T,
          repel = T, label.size = 2) +
  ggtitle("DC & macrophage in epidermis") +
  theme_jpub +
  NoLegend()

publish_pdf('mission/SLE_TRPM2_MfMo/figs/epidermis.dc.mf.umap.pdf')

mono.lc.marker <- c('LYZ','IL1B','CLEC10A','CD1C',
                    'S100B','CD1A','CD207','FCGBP')

sobj.dc |> DotPlot(mono.lc.marker, cols = 'RdYlBu')

visual.mo <- c(2,4,8,9,10,13,14)
visual.lc <- c(1,3,5,7,11)

sobj.dc <- sobj.dc |>
  mutate(manual.fine = case_when(seurat_clusters %in% visual.lc ~ 'LC',
                                 seurat_clusters %in% visual.mo ~ 'MF',
                                 .default = 'other'))

sobj.dc |> DimPlot(group.by = 'manual.fine')

sobj.dc <- sobj.dc |>
  mutate(m2.fine = str_c(manual.fine, seurat_clusters, sep = '_'))

sobj.dc |> FeaturePlot('TRPM2', split.by = 'group', order = T,
                       cols = c('lightgrey','red'))

### M2 expr ------
sobj.dc |>
  filter(group == 'SLE') |>
  DotPlot('TRPM2', cols = 'RdYlBu', group.by = 'm2.fine')

lcmo.m2.expr <- last_plot() |> pluck('data')

lcmo.m2.expr |> ggplot(aes(id, avg.exp)) + geom_point()

lcmo.m2.expr |>
  filter(str_detect(id, 'MF')) |>
  mutate(kmeans)

### frac SLEvHC -------
lcmo.leiden.frac <- sobj.dc |>
  calc_frac_conf_on_grouped_count(group, m2.fine)

lcmo.leiden.frac |>
  ggplot(aes(group, m2.fine, size = fraction)) +
  geom_point() +
  geom_label(aes(label = n), nudge_y = .3, size = 4)

lcmo.leiden.delta <- lcmo.leiden.frac |>
  pivot_wider(names_from = group, values_from = fraction,
              values_fill = 0, id_cols = m2.fine) |>
  mutate(delta = SLE - HC)

lcmo.leiden.delta |>
  left_join(lcmo.m2.expr, join_by(m2.fine == id)) |> DT::datatable()
  ggplot(aes(avg.exp, delta)) +
  geom_smooth(method = 'lm') +
  geom_point() +
  stat_cor()

### M2 SLEvHC ---------
lcmo.m2.slevhc <- sobj.dc |>
  FindMarkersAcrossVar(split.by = 'm2.fine', group.by = 'group',
                       ident.1 = 'SLE', features = 'TRPM2')

lcmo.m2.slevhc

m2h.excel <- c(12,4,10,2)

sobj.dc <- sobj.dc |>
  mutate(m2.main = case_when(str_detect(m2.fine, 'LC') ~ 'DC',
                             seurat_clusters %in% m2h.excel ~ 'M2-hi MF',
                             seurat_clusters %in% c(9,13) ~ 'M2-mid MF',
                             .default = 'M2-lo MF'))

#### M2 dotplot in epi MF -------
sobj.dc |>
  filter(m2.main != 'LC', group == 'SLE', seurat_clusters != 12) |>
  DotPlot('TRPM2') +
  scale_color_gradient2(low = kn.pal2[2], high = kn.pal2[1]) +
  labs(x = 'Gene', y = 'Macrophage clusters')

publish_pdf('mission/SLE_TRPM2_MfMo/figs/sle.epidermis.mf.m2.dotplot.pdf',
            width = 80, height = 75)

last_plot() |>
  pluck('data') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/sle.epidermis.mf.m2.dotplot.csv')

sobj.dc |> DimPlot(group.by = 'm2.main')

#### frac change ----------
sobj.dc |>
  calc_frac_conf_on_grouped_count(orig.ident, m2.main) |>
  mutate(group = ifelse(str_detect(orig.ident, 'SLE'), 'SLE', 'HC')) |>
  ggplot(aes(group, fraction)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~m2.main, scales = 'free_y') +
  stat_compare_means(method = 't.test', comparisons = list(c('SLE','HC'))) +
  labs(y = '% in DC & macrophages', title = 'Macrophage & DC in SLE skin') +
  theme_bw()

last_plot() |>
  pluck('data') |>
  select(m2.main, fraction, group) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/sle.epidermis.MF.3subsets.percentage.csv')

sobj <- 
sobj.dc@meta.data |>
  as_tibble(rownames = '.cell') |>
  select(.cell, m2.main) |>
  left_join(sobj, y = _)

sobj |>
  calc_frac_conf_on_grouped_count(orig.ident, m2.main) |>
  drop_na() |>
  mutate(group = ifelse(str_detect(orig.ident, 'SLE'), 'SLE', 'HC')) |>
  ggplot(aes(group, fraction, color = group)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~m2.main, scales = 'free_y') +
  stat_compare_means(method = 't.test', comparisons = list(c('SLE','HC'))) +
  labs(y = '% in epidermis cells', title = 'Macrophage & DC in SLE epidermis') +
  theme_bw() +
  scale_color_manual(values = kn.pal2)

sobj.dc |>
  FindMarkersAcrossVar(split.by = 'm2.main', group.by = 'group',
                       ident.1 = 'SLE', features = 'TRPM2',
                       min.pct = 0, logfc.threshold = 0)

## M2 expression in all cell types ------
sobj |>
  mutate(sle_tissue = str_c(tissue, '_', sle)) |>
  DotPlot2d('TRPM2', sle_tissue, manual_main) +
  labs(title = 'TRPM2 expression in SLE skin tissue',
       x = 'Sample', y = 'Cell type')

g1 <- last_plot()

g1$data |>
  mutate(group.y = fct_reorder(group.y, avg.exp)) |>
  BubblePlot(d2 = T) +
  labs(x = 'Group', y = 'Cell type', title = 'TRPM2 expression in skin') +
  theme_jpub +
  RotatedAxis()

publish_source_plot('skin.celltype.trpm2.dotplot', project = proj.nm,
                    height = 55, width = 60)

## M2 expr in macro by sample ----------
mf.m2.sle.ts.samp <- sobj |>
  filter(manual_main == 'Macrophage') |>
  AverageExpression(features = 'TRPM2',
                    group.by = c('sle', 'tissue', 'orig.ident'))

mf.m2.sle.ts.samp$RNA |>
  as_tibble() |>
  pivot_longer(everything()) |>
  add_row(name = rep.int('HC_epidermis_1', 4), value = rep.int(0, 4)) |>
  mutate(tissue = str_extract(name, '(epi|)dermis'),
         group = str_extract(name, 'HC|SLE')) |>
  tidyplot(x = tissue, y = value, color = group) |>
  add_boxplot() |>
  add_data_points_beeswarm() |>
  add_test_pvalue(hide_info = T) +
  labs(title = 'TRPM2 expression', y = 'Normalized expression') 

save_plot(plot = last_plot(),'sle.hc.epi.derm.trpm2.expr.pdf')

## cell fraction ----------
type.conf <- sobj |>
  calc_frac_conf_on_grouped_count(sle, manual_main)

type.conf |>
  filter(manual_main %in% c('Macrophage','DC','B_cell')) |>
  ggplot(aes(sle, fraction)) +
  geom_col() +
  facet_wrap(~manual_main)

sample.conf <- sobj |>
  calc_frac_conf_on_grouped_count(orig.ident, manual_main)

sample.conf <- sobj |>
  as_tibble() |>
  distinct(orig.ident, sle, tissue) |>
  right_join(sample.conf) 

sample.conf |>
  filter(manual_main %in% c('Macrophage','DC','B_cell')) |>
  ggplot(aes(x = sle, y = fraction * 100, fill = sle, group = sle)) +
  stat_summary(geom = 'col', fun = mean) +
  geom_dotgraph +
  facet_wrap(vars(tissue, manual_main), scales = 'free_y') +
  theme_bw() +
  labs(x = 'Group', y = '% in tissue', fill = 'Group',
       title = 'Immune cells in SLE dermis and epidermis') +
  scale_fill_hue(direction = -1)

## epidermis ---------
sobj.epi <- read_rds('mission/FPP/xiangya_sle_scRNA/mix4.epi.sle.rds')

sobj |> DimPlot(group.by = 'manual.main')

sobj |> DimPlot(group.by = 'seurat_clusters')

sobj.myl <- sobj |> filter(manual.main == 'DC')

sobj.myl <- sobj.myl |>
  quick_process_seurat(c('orig.ident', 'batch'), skip_norm = T)

sobj.myl |> DotPlot('TRPM2', cols = 'RdYlBu')

sobj.myl |> FindAllMarkers(features = 'TRPM2')

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj.myl <- sobj.myl |>
  mark_cell_type_singler(hpca, new_label = 'hpca.main')

sobj.myl |>
  DotPlot2d('TRPM2', group, hpca.main)

sobj.myl |>
  dplyr::count(hpca.main, group)

sobj.myl |>
  filter(hpca.main != 'B_cell') |>
  FindMarkersAcrossVar(split.by = 'hpca.main', group.by = 'group',
                       ident.1 = 'SLE', features = 'TRPM2')

epi.mf.m2 <- sobj.myl |>
  filter(hpca.main == 'Monocyte') |>
  AverageExpression(features = 'TRPM2', group.by = 'orig.ident')

epi.mf.m2$RNA |>
  as_tibble() |>
  pivot_longer(everything()) |>
  mutate(tissue = ifelse(str_detect(name, 'D$'), 'Dermis', 'Epidermis'),
         group = str_extract(name, 'SLE|HC'))
sobj |>
  DotPlot2d('TRPM2', group, manual.fine) +
  labs(title = 'TRPM2 expression in SLE skin tissue',
       x = 'Sample', y = 'Cell type')

sobj |>
  DotPlot(seurat_markers, cols = 'RdYlBu')

g1 <- last_plot()

g1$data |>
  mutate(group.y = fct_reorder(group.y, avg.exp)) |>
  BubblePlot(d2 = T) +
  labs(x = 'Group', y = 'Cell type', title = 'TRPM2 expression in skin')

## M2-hi subsets in macrophage? -------
sobj.mf <- sobj |> filter(manual_main == 'Macrophage')

sobj.mf <- sobj.mf |>
  quick_process_seurat(c('orig.ident', 'dataset', 'tissue'), skip_norm = T)

### SAVE rds --------
sobj.mf |> write_rds('mission/SLE_TRPM2_MfMo/data/sle.skin.mf.rds')

sobj.mf <- read_rds('mission/SLE_TRPM2_MfMo/data/sle.skin.mf.rds')

### dotplot --------
sobj.mf |> DotPlot(c('TRPM2','CD14','FCGR3A'), cols = 'RdBu') +
  labs(x = 'Gene', y = 'Macrophage clusters')

publish_pdf('mission/SLE_TRPM2_MfMo/figs/skin.mf.m2.dotplot.pdf',
            width = 85, height = 88)

sobj.mf$m2.type |> unique()

m2h_bc <- sobj.mf |>
  filter(m2.type == 'TRPM2-high MF') |>
  colnames()

sobj <- sobj |>
  mutate(manual_fine = case_when(.cell %in% m2h_bc ~ 'CD14+ Macrophage',
                                 manual_main == 'Macrophage' ~ 'CD14- Macrophage',
                                 .default = manual_main))

sobj |>
  DotPlot2d('TRPM2', sle, manual_fine)

g1 <- last_plot()

g1$data |>
  mutate(group.y = fct_reorder(group.y, avg.exp, max)) |>
  BubblePlot(d2 = T) +
  labs(x = 'Group', y = 'Cell type', title = 'SLE skin TRPM2') +
  theme_jpub()

publish_source_plot('sle.skin.cd14.mf.m2.dotplot', width = 60)

### module score of inflammatory genes ----------
flare.gsebp.tb <-
  read_csv('mission/SLE_TRPM2_MfMo/flare.m2hvl.gsebp.csv')

inflam.gene <- flare.gsebp.tb |>
  filter(str_detect(Description, '^inflamma')) |>
  pull(core_enrichment) |>
  str_split_1('/')

sobj.mf <- sobj.mf |>
  AddModuleScore(features = list(inflam.gene), name = 'inflam')

mf.infl <- sobj.mf |>
  summarise(inflam.score = mean(inflam1), .by = seurat_clusters)

### define M2-hi MF --------------
sobj.mf |>
  filter(sle == 'HC') |>
  FindAllMarkers(features = 'TRPM2')

sobj.mf |> DimPlot(cols = 'Paired', label = T, label.box = T,
                   label.size = 2, repel = T) +
  theme_jpub

publish_pdf('mission/SLE_TRPM2_MfMo/figs/skin.mf.umap.pdf', 60, 50)

sobj.mf |>
  filter(seurat_clusters != 2) |>
  DotPlot2d('TRPM2', sle, seurat_clusters)

sobj.mf |>
  DotPlot(c('TRPM2','CD14','CCR2'), cols = c('RdYlBu'))

mf.m2.exp <- last_plot() |>
  pluck('data')

mf.m2.exp |>
  as_tibble() |>
  mutate(subtype = case_when(id %in% c(3,8) ~ 'TRPM2-high',
                             id %in% c(4,6) ~ 'TRPM2-mid',
                             .default = 'TRPM2-low') |>
           fct_relevel('TRPM2-high', 'TRPM2-mid')) |>
  left_join(mf.infl, join_by(id == seurat_clusters)) |>
  ggplot(aes(avg.exp, inflam.score)) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2) +
  theme_bw() +
  stat_cor(size = 2) +
  labs(x = 'Average expression of TRPM2', y = 'Inflammatory signature score',
       title = 'SLE skin macrophage clusters') +
  theme_jpub

mf.m2.exp |>
  as_tibble() |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  mutate(subtype = case_when(id %in% c(3,8) ~ 'TRPM2-high',
                             id %in% c(4,6) ~ 'TRPM2-mid',
                             .default = 'TRPM2-low') |>
           fct_relevel('TRPM2-high', 'TRPM2-mid')) |>
  ggplot(aes(CCR2, CD14)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2) +
  theme_bw() +
  stat_cor(size = 2) +
  labs(title = 'SLE skin macrophage clusters') +
  theme_jpub()

publish_source_plot('sle.skin.mf.ccr2.cd14.score.correlation', width = 70)

sle_skin_mf_tag <-
  read_csv('mission/SLE_TRPM2_MfMo/results/sle.skin.mf.ccr2.cd14.score.correlation.csv')

sle_skin_mf_tag |>
  ggplot(aes(TRPM2, CCR2)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2) +
  stat_cor(size = 2) +
  labs(title = 'SLE skin macrophage') +
  theme_jpub()

publish_source_plot('sle.skin.mf.ccr2.trpm2.score.correlation', width = 70)

sle.skin.mf.m2.inflam.score.correlation |>
  ggplot(aes(avg.exp, inflam.score)) +
  geom_smooth(method = 'lm', linetype = 'dashed', color = 'grey') +
  geom_point(aes(color = subtype)) +
  geom_text_repel(aes(label = id), size = 2) +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  stat_cor(size = 2) +
  labs(x = 'Average expression of TRPM2', y = 'Inflammatory signature score',
       title = 'SLE skin macrophage clusters') +
  theme_jpub

publish_source_plot('sle.skin.mf.m2.inflam.score.correlation', width = 70)

sle.skin.mf.m2.inflam.score.correlation <-
  read_csv('mission/SLE_TRPM2_MfMo/results/sle.skin.mf.m2.inflam.score.correlation.csv')

mf.m2.exp |>
  ggplot(aes(avg.exp)) +
  geom_density()

mf.m2.exp |>
  mutate(cluster = kmeans(avg.exp, 2)$cluster) |>
  as_tibble()

sobj.mf <- sobj.mf |>
  mutate(m2.type = ifelse(seurat_clusters %in% c(3,8), 'TRPM2-high MF',
                          'TRPM2-low MF'))

sobj.mf |> DimPlot(split.by = 'tissue', group.by = 'm2.type') +
  ggtitle('Macrophages in SLE skin')

#### cell frac logfc ----------
#### in macrophage --------
mf.leiden.frac <- sobj.mf@meta.data |>
  as_tibble(rownames = '.cell') |>
  discov_frac_change(sle, seurat_clusters, SLE, HC)

mf.m2.exp |>
  as_tibble() |>
  left_join(mf.leiden.frac, join_by(id == subtype)) |>
  ggplot(aes(avg.exp, log2fc_frac)) +
  geom_hline(yintercept = 0) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = type)) +
  geom_text_repel(aes(label = id)) +
  theme_bw() +
  stat_cor(label.x.npc = 'middle') +
  scale_color_manual(values = c('blue','red','green')) +
  labs(x = 'Average expression of TRPM2', y = 'Log2FC of % in MF SLE/HC',
       title = 'SLE skin macrophage clusters', color = '% in MF SLE/HC')

#### in (epi)dermis --------
mf.dermi.leiden.frac <- sobj.mf |>
  filter(tissue == 'dermis') |>
  discov_frac_change(sle, seurat_clusters, SLE, HC)

mf.m2.exp |>
  as_tibble() |>
  left_join(mf.dermi.leiden.frac, join_by(id == subtype)) |>
  ggplot(aes(avg.exp, log2fc_frac)) +
  geom_hline(yintercept = 0) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = type)) +
  geom_text_repel(aes(label = id)) +
  theme_bw() +
  stat_cor(label.x.npc = 'middle') +
  scale_color_manual(values = c('blue','red','green')) +
  labs(x = 'Average expression of TRPM2', y = 'Log2FC of % in dermal MF SLE/HC',
       title = 'SLE skin macrophage clusters', color = '% in MF SLE/HC')

mf.epi.leiden.frac <- sobj.mf |>
  filter(tissue == 'epidermis') |>
  discov_frac_change(sle, seurat_clusters, SLE, HC)

mf.m2.exp |>
  as_tibble() |>
  left_join(mf.epi.leiden.frac, join_by(id == subtype)) |>
  ggplot(aes(avg.exp, log2fc_frac)) +
  geom_hline(yintercept = 0) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = type)) +
  geom_text_repel(aes(label = id)) +
  theme_bw() +
  stat_cor(label.x.npc = 'middle') +
  scale_color_manual(values = c('blue','red','green')) +
  labs(x = 'Average expression of TRPM2',
       y = 'Log2FC of % in epidermal MF SLE/HC',
       title = 'SLE skin macrophage clusters', color = '% in MF SLE/HC')

#### in total skin ----------
mf.leiden <- sobj.mf@meta.data |>
  as_tibble(rownames = '.cell') |>
  mutate(.cell, mf.leiden = seurat_clusters, .keep = 'none')

skin.leiden.frac <- sobj@meta.data |>
  as_tibble(rownames = '.cell') |>
  left_join(mf.leiden) |>
  discov_frac_change(sle, mf.leiden, SLE, HC)

mf.m2.exp |>
  as_tibble() |>
  left_join(skin.leiden.frac, join_by(id == subtype)) |>
  ggplot(aes(avg.exp, log2fc_frac)) +
  geom_hline(yintercept = 0) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = type)) +
  geom_text_repel(aes(label = id)) +
  theme_bw() +
  stat_cor(label.x.npc = 'middle') +
  scale_color_manual(values = c('blue','red','green')) +
  labs(x = 'Average expression of TRPM2', y = 'Log2FC of % in skin SLE/HC',
       title = 'SLE skin macrophage clusters', color = '% in skin SLE/HC')

sobj |>
  calc_frac_conf_on_grouped_count(sle, manual_main)

#### M2 expr SLE vs HC --------
mf.leiden.m2exp.slevhc <- sobj.mf |>
  FindMarkersAcrossVar(split.by = 'seurat_clusters', group.by = 'sle',
                       ident.1 = 'SLE', features = 'TRPM2', min.pct = 0,
                       logfc.threshold = 0)

mf.m2.exp |>
  as_tibble() |>
  inner_join(mf.leiden.m2exp.slevhc, join_by(id == cluster)) |>
  ggplot(aes(avg.exp, avg_log2FC)) +
  geom_hline(yintercept = 0) +
  geom_smooth(method = 'lm', linetype = 'dashed', se = FALSE, color = 'grey') +
  geom_point(aes(color = avg_log2FC > 0)) +
  geom_text_repel(aes(label = id)) +
  theme_bw() +
  stat_cor(label.x.npc = .3) +
  scale_color_manual(values = c('blue','red'),
                     labels = c('Downregulated','Upregulated')) +
  labs(x = 'Average expression of TRPM2', y = 'Log2FC of TRPM2 SLE/HC',
       title = 'SLE skin macrophage clusters', color = 'TRPM2 in SLE/HC')

#### TRPM2-mid MF ------------
sobj.mf$m2.type |> table()

sobj.mf <- sobj.mf |>
  mutate(m2.mid = ifelse(seurat_clusters %in% c(4,6), 'TRPM2-mid MF', m2.type) |> fct_relevel('TRPM2-high MF', 'TRPM2-mid MF'))

sobj.mf$m2.mid |> table()

### TRPM family -----------
trpmx <- sobj.mf |> rownames() |> str_subset('TRPM\\d') |> str_sort()

trpmx.slevhc <- sobj.mf |>
  FindMarkers(group.by = 'sle', ident.1 = 'SLE',
              features = trpmx, min.pct = 0, logfc.threshold = 0)

trpmx.slevhc |>
  as_tibble(rownames = 'gene') |>
  mutate(logp = -log10(p_val) * avg_log2FC / abs(avg_log2FC)) |>
  ggplot(aes(logp, gene, fill = avg_log2FC)) +
  geom_col() +
  geom_vline(xintercept = 0) +
  geom_vline(xintercept = c(-log10(.05),log10(.05)), linetype = 'dashed') +
  scale_fill_gradient2(low = 'blue', high = 'red',mid = '#ddd') +
  theme_pubr(legend = 'right') +
  labs(x = '-log10(adjusted P value)',
       title = 'TRPM family in macrophages from SLE skin',
       subtitle = 'GSE179633 (n=17)')

### compare cell frac ----------
skin.m2hl.conf <- sobj |>
  calc_frac_conf_on_grouped_count(orig.ident, m2.type)

sobj.epi <- read_rds('mission/FPP/xiangya_sle_scRNA/mix4.epi.sle.rds')

sobj.epidc <- sobj.epi |>
  filter(manual.main == 'DC') |>
  quick_process_seurat(skip_norm = T)

sobj.epidc |>
  DotPlot('TRPM2', cols = 'RdYlBu')

m2.epidc <- sobj.epidc |>
  filter(seurat_clusters %in% c(6,8)) |>
  pull(.cell)

epin.conf <- sobj.epi |>
  filter(sobj.epi$batch == 'nakamizu') |>
  mutate(m2.type = .cell %in% m2.epidc) |>
  calc_frac_conf_on_grouped_count(orig.ident, m2.type) |>
  filter(m2.type) |>
  slice_min(fraction, n=4) |>
  mutate(tissue = 'Epidermis', group = 'HC', m2.type = 'm2hi',
         fraction = fraction/1)

m2hl.conf |>
  mutate(tissue = ifelse(str_detect(orig.ident, 'D$'),
                         'Dermis', 'Epidermis'),
         group = str_extract(orig.ident, 'SLE|HC')) |>
  filter(str_detect(m2.type, 'high')) |>
  bind_rows(epin.conf) |>
  tidyplot(x = tissue, y = fraction, color = group) |>
  add_sem_errorbar() |>
  add_mean_bar() |>
  add_data_points_beeswarm(color = 'black', show.legend = F) |>
  add_test_pvalue(hide_info = T) +
  labs(x = 'Tissue', y = 'Fraction in skin cells',
       title = 'Proportion of TRPM2-high macrophages')

last_plot() |>
  save_plot('mission/SLE_TRPM2_MfMo/figs/SLE.hc.skin.m2.frac.pdf')

### TRPM2 expr in SLE/HC & hi/mid/lo ------
sobj.mf |>
  bill.violin('TRPM2', m2.mid)

m2.percell <- sobj.mf |>
  get_abundance_sc_wide('TRPM2')

sobj.mf |>
  left_join(m2.percell) |>
  ggplot(aes(sle, TRPM2)) +
  geom_violin(aes(fill = sle), scale = 'width') +
  stat_summary(fun = logtpm.mean, geom = 'crossbar',
               width = .3, color = 'black') +
  facet_wrap(~m2.mid) +
  theme_pubr()

mf.m2mid.slehc.m2 <- last_plot() |>
  pluck('data')

mf.m2mid.slehc.m2 |>
  select(.cell, sle, m2.mid, TRPM2) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/skin.m2.mf.3subsets.sle-hc.m2.violin.csv')

mf.m2mid.slehc.m2 <-
  read_csv('mission/SLE_TRPM2_MfMo/results/skin.m2.mf.3subsets.sle-hc.m2.violin.csv')

mf.m2mid.slehc.m2 |>
  mutate(m2.mid = fct_relevel(m2.mid, 'TRPM2-high MF', 'TRPM2-mid MF')) |>
  ggplot(aes(sle, TRPM2)) +
  geom_violin(aes(fill = sle), scale = 'width') +
  stat_summary(fun = logtpm.mean, geom = 'crossbar',
               width = .3, color = 'black') +
  facet_wrap(~m2.mid) +
  theme_pubr() +
  labs(x = 'Group', y = 'Normalized expression', fill = 'Group',
       title = 'TRPM2 expression in skin macrophages') +
  scale_fill_hue(direction = -1)

sobj.mf |>
  FindMarkersAcrossVar(split.by = 'm2.mid', group.by = 'sle',
                       ident.1 = 'SLE', features = 'TRPM2')

### TRPM2 expr hi vs lo ---------
mf.m2.expr <- sobj.mf |>
  bill.violin('TRPM2', m2.type) |>
  pluck('data')

mf.m2.expr |>
  tidyplot(x = m2.type, y = .abundance_RNA, color = m2.type) |>
  add_violin() |>
  add_mean_dash() |>
  adjust_legend_title('') |>
  reorder_x_axis_labels('TRPM2-low MF') |>
  add_test_pvalue(method = 'wilcox.test', hide_info = T) |>
  adjust_y_axis(limits = c(0,3.3)) +
  labs(x = 'Cell type', y = 'Normalized expression',
       title = 'TRPM2 expression in macrophage')

last_plot() |>
  save_plot('SLE.skin.m2hvl.m2.expr.pdf')

sobj.mf |>
  FindMarkersAcrossVar(split.by = 'seurat_clusters', group.by = 'sle',
                       ident.1 = 'SLE', features = 'TRPM2')

m2.expr <- sobj.mf |>
  AverageExpression(features = 'TRPM2', group.by = c('orig.ident'))

m2.expr <- m2.expr$RNA |>
  as_tibble() |>
  pivot_longer(everything()) |>
  mutate(tissue = ifelse(str_detect(name, 'D$'), 'Dermis', 'Epidermis'),
         group = str_extract(name, 'SLE|HC'))

m2.expr |>
  filter(tissue == 'Dermis') |>
  ggplot(aes(x = group, y = value, color = group)) +
  geom_boxplot() +
  geom_dotplot(binaxis = 'y', stackdir = 'center', fill = 'white') +
  labs(title = 'TRPM2 in macrophage from dermis', y = 'Normalized expression') +
  theme_pubr() +
  scale_color_hue(direction = -1)

### DEGA m2h vs m2l -----
library(clusterProfiler)
library(enrichplot)
m2hvl.slehc.deg <- sobj.mf |>
  FindMarkersAcrossVar(split.by = 'sle', group.by = 'm2.type',
                       ident.1 = 'TRPM2-high MF')

m2hvl.slehc.deg |>
  write_csv('mission/SLE_TRPM2_MfMo/results/m2hvl.skin.deg.csv')

m2hvl.slehc.deg <-
  read_csv('mission/SLE_TRPM2_MfMo/results/m2hvl.skin.deg.csv')

sle.m2hvl.gse <- m2hvl.slehc.deg |>
  filter(cluster == 'SLE', p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL')

sle.m2hvl.gse@result |>
  as_tibble() |>
  mutate(tags = str_extract(leading_edge, '(?<=tags=)\\d+') |> as.numeric(),
         signal = str_extract(leading_edge, '(?<=signal=)\\d+') |> as.numeric()) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/sle.skin.m2hvl.gsego.csv')

sle.m2hvl.gse.res <-
  read_csv('mission/SLE_TRPM2_MfMo/results/sle.skin.m2hvl.gsego.csv')

set.seed(42)

sle.m2hvl.gse.res |>
  filter(NES > 0) |>
  plot_enrichment(force_regex = 'inf') +
  theme_jpub() +
  labs(x = 'Normalized Enrichment Score', y = 'Pathway',
       title = 'GO BP pathway enrichment in SLE skin:\nTRPM2-hi macro vs TRPM2-mid/lo macro')

publish_source_plot('SLE.skin.m2h.vs.m2ml.mf.gsego', width = 70)

sle.m2hvl.gse@result |>
  filter(NES > 0, qvalue < .05, ONTOLOGY == 'BP') |>
  DT::datatable()

sle.m2hvl.gse |>
  gseaplot2('GO:0050729', base_size = 8,
            title = 'GO pathway: Positive regulation of inflammatory response')

sle.m2hvl.gse <- sle.m2hvl.gse |>
  simplify()

infla.path <- sle.m2hvl.gse@result |>
  filter(ID == 'GO:0050729')

sle.m2hvl.gse@result |>
  filter(NES > 1, qvalue < .05, ONTOLOGY == 'BP') |>
  head(10) |>
  mutate(Description = str_wrap(Description, 50)) |>
  tidyplot(x = NES, y = Description, color = p.adjust) |>
  add_mean_bar() |>
  sort_y_axis_labels() +
  scale_fill_distiller(palette = 'Reds') +
  labs(x = 'Normalized Enrichment Score', y = 'Pathway',
       title = 'GO BP pathway enrichment in SLE skin:\nTRPM2-hi macro vs TRPM2-mid/lo macro') 

last_plot() |>
  save_plot('mission/SLE_TRPM2_MfMo/figs/SLE.skin.m2hvl.gogse.pdf')

sle.m2hvl.gse@result |>
  filter(NES > 1, qvalue < .05, ONTOLOGY == 'BP') |>
  slice_sample(n = 9) |>
  bind_rows(infla.path) |>
  distinct(ID, .keep_all = T) |>
  mutate(Description = str_wrap(Description, 50)) |>
  tidyplot(x = NES, y = Description, color = p.adjust) |>
  add_mean_bar() |>
  sort_y_axis_labels() +
  scale_fill_distiller(palette = 'Reds') +
  labs(x = 'Normalized Enrichment Score', y = 'Pathway',
       title = 'GO BP pathway enrichment in SLE skin:\nTRPM2-hi macro vs TRPM2-mid/lo macro')

last_plot() |> 
  save_plot('mission/SLE_TRPM2_MfMo/figs/SLE.skin.m2hvl.gogse.pdf')

### Top chemotaxis genes --------
skin.immune <- sle.m2hvl.gse@result |>
  filter(Description == 'positive regulation of immune system process') |>
  pull(core_enrichment) |>
  str_split_1('/')

m2hvl.slehc.deg |>
  filter(cluster == 'SLE', gene %in% skin.immune) |>
  slice_min(p_val_adj, n = 20) |>
  mutate(`-log10(p_val_adj)` = -log10(p_val_adj)) |>
  tidyplot(x = avg_log2FC, y = gene, color = `-log10(p_val_adj)`) |>
  add_mean_bar() |>
  sort_y_axis_labels() +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  ggtitle('Top enriched immune response genes in SLE skin TRPM2-hi vs TRPM2-lo')

last_plot() |> 
  save_plot('SLE.skin.m2hvl.immune.genes.pdf')

m2.top15 <- m2hvl.slehc.deg |>
  filter(cluster == 'SLE', gene %in% skin.immune) |>
  slice_min(p_val_adj, n = 15) |>
  pull(gene)

sobj.mf |>
  filter(sle == 'SLE') |>
  bill.violin(m2.top15, m2.mid, facet.ncol = 5) +
  labs(x = 'Cell type', y = 'Normalized expression', fill = 'Cell type',
       title = '"Positive regulation of immune system" genes in SLE skin') +
  theme_jpub +
  theme(axis.text.x = element_blank())

publish_pdf('sle.skin.mf.3subsets.positive.immune.violin.pdf', width = 90)

## top inflamm genes --------
top15.inflam <-
  read_csv('mission/SLE_TRPM2_MfMo/results/top15.inflamm.gene.csv')

sle.m2hvl.gse.res <-
  read_csv('mission/SLE_TRPM2_MfMo/results/sle.skin.m2hvl.gsego.csv')

sle.infl.gene <- sle.m2hvl.gse.res |>
  filter(str_detect(Description, 'infl')) |>
  pull(core_enrichment) |>
  str_split_1('/')

skin.m2h.infl14 <- m2hvl.slehc.deg |>
  filter(gene %in% sle.infl.gene, cluster == 'SLE') |>
  slice_min(p_val_adj, n = 14)

sobj.mf |>
  filter(sle == 'SLE') |>
  bill.violin(c('CD14', skin.m2h.infl14$gene),
              m2.type, facet.ncol = 5, pubsize = T) +
  labs(x = 'Cell type', y = 'Normalized expression', fill = 'Cell type',
       title = 'Inflammatory response genes in SLE skin')

publish_pdf('sle.skin.mf.2subsets.inflamm.violin.pdf', width = 90)

## heatmap of m2-hi/lo & DC -----------
m2hvl.deg <- sobj.mf |>
  FindMarkers(group.by = 'm2.type', ident.1 = 'TRPM2-low MF')

m2l.mrk <- m2hvl.deg |>
  as_tibble(rownames = 'gene') |>
  filter(p_val_adj < .01, avg_log2FC > 1) |>
  slice_min(p_val_adj, n = 30) |>
  pull(gene)

m2h.mrk <- m2hvl.deg |>
  as_tibble(rownames = 'gene') |>
  filter(p_val_adj < .01, avg_log2FC < -1) |>
  slice_min(p_val_adj, n = 30) |>
  pull(gene)
  
g1 <- sobj.mf |>
  ScaleData(features = c(m2h.mrk, m2l.mrk)) |>
  DoHeatmap(c(m2h.mrk, m2l.mrk), group.by = 'm2.mid', label = F,
            raster = F)

g1[[1]] +
  theme(axis.text.y = element_text(size = 2),
        text = element_text(size = 5))

publish_pdf('sle.skin.mf.2subsets.heatmap.pdf', width = 90)

sobj.mf |>
  ScaleData(features = c(m2h.mrk, m2l.mrk)) |>
  DoHeatmap(c(m2h.mrk, m2l.mrk), group.by = 'seurat_clusters')

## m2hi/lo DC ---------
sobj.m2dc <- sobj |>
  filter(manual_main == 'DC')

sobj.m2dc <- sobj.m2dc |>
  quick_process_seurat(batch = c('orig.ident','dataset'), skip_norm = T)

### UMAP ---------
sobj.m2dc |>
  DimPlot(cols = paired.pal, label = T, repel = T,
          label.size = 1.5, label.box = T) +
  ggtitle('Skin DC') +
  theme_jpub +
  NoLegend()

publish_pdf('mission/SLE_TRPM2_MfMo/figs/skin.dc.umap.pdf',
            50, height = 50)

### M2 dotplot --------
sobj.m2dc |>
  DotPlot('TRPM2') +
  scale_color_gradient2(low = kn.pal2[2], high = kn.pal2[1]) +
  labs(x = 'Gene', y = 'DC clusters', title = 'Skin DC')

publish_pdf('mission/SLE_TRPM2_MfMo/figs/skin.dc.trpm2.dotplot.pdf',
            85, height = 90)

sobj.m2dc <- sobj.m2dc |>
  get_abundance_sc_wide('TRPM2') |>
  left_join(x = sobj.m2dc, y = _)

### define m2h DC --------
sobj.m2dc <- sobj.m2dc |>
  mutate(m2dc.main = ifelse(seurat_clusters %in% c(2,4,7,14),
                            'TRPM2-high DC', 'TRPM2-low DC'))

### frac change -------
m2dc.frac <- sobj.m2dc@meta.data |>
  as_tibble(rownames = '.cell') |>
  select(.cell, m2dc.main) |>
  left_join(x = sobj, y = _) |>
  calc_frac_conf_on_grouped_count(orig.ident, m2dc.main)

m2dc.frac |>
  mutate(group = str_extract(orig.ident, 'HC|SLE')) |>
  filter(!is.na(m2dc.main)) |>
  select(group, fraction, m2dc.main) |>
  ggplot(aes(group, fraction)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~m2dc.main) +
  stat_compare_means(method = 't.test',
                     comparisons = list(c('SLE','HC'))) +
  theme_jpub

publish_source_plot('sle.skin.dc.subset.frac', proj.nm)

### M2 violin -----------
sobj.m2dc |>
  select(sle, TRPM2, m2dc.main) |>
  ggplot(aes(sle, TRPM2)) +
  geom_violin(aes(fill = sle), scale = 'width') +
  stat_summary(fun = logtpm.mean, geom = 'crossbar',
               width = .3, color = 'black') +
  facet_wrap(~m2dc.main) +
  theme_pubr() +
  labs(x = 'Group', y = 'Normalized expression', fill = 'Group',
       title = 'TRPM2 expression in skin DC') +
  scale_fill_hue(direction = -1)

publish_source_plot('skin.dc.subset.slevhc.m2.violin', project = proj.nm,
                    80, 80)

m2dc.slevhc.deg <- sobj.m2dc |>
  FindMarkersAcrossVar(split.by = 'm2dc.main', group.by = 'sle',
                       ident.1 = 'SLE')

m2dc.slevhc.deg |>
  filter(gene == 'TRPM2')

m2hdc.slevhc.lst <- m2dc.slevhc.deg |>
  filter(cluster == 'TRPM2-high DC', p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T)

m2hdc.gsego <- m2hdc.slevhc.lst |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL')

m2hdc.gsego@result |>
  filter(qvalue < .05, str_length(Description) < 60) |>
  slice_sample(n = 10) |>
  mutate(Description = str_wrap(Description, 50)) |>
  tidyplot(x = NES, y = Description, color = p.adjust) |>
  add_mean_bar() |>
  sort_y_axis_labels() +
  scale_fill_distiller(palette = 'Reds') +
  labs(x = 'Normalized Enrichment Score', y = 'Pathway',
       title = 'GO BP pathway enrichment in SLE skin:\nTRPM2-hi DC vs TRPM2-mid/lo DC')

publish_source_plot('sle.skin.m2h.dc.go.gsea', project = proj.nm,
                    width = 100, height = 60)

sobj.m2dc |>
  bill.violin(top15.inflam$gene, group.by = m2dc.main, facet.ncol = 5) +
  labs(x = 'Cell type', y = 'Normalized expression', fill = 'Cell type',
       title = 'Inflammatory response genes in SLE skin DC') +
  theme_jpub +
  theme(axis.text.x = element_blank())

publish_pdf('sle.skin.dc.inflamm.violin.pdf', project = proj.nm,
            width = 90)

m2dc.inflam <- last_plot() |>
  pluck('data')

logmean <- m2dc.inflam |>
  summarise(lmean = logtpm.mean(.abundance_RNA), .by = c(m2dc.main, .feature)) |>
  rename('Group' = 'm2dc.main')

m2dc.inflam |>
  rename('Normalized expression' = '.abundance_RNA',
         'Group' = 'm2dc.main') |>
  tidyplot(Group, `Normalized expression`, color = Group) |>
  add_violin(trim = T) |>
  add_mean_dash(data = logmean, aes(y = lmean)) |>
  remove_x_axis_labels() |>
  add_title('Inflammatory response genes in SLE skin DC') |>
  adjust_colors(new_colors = kn.pal2) |>
  split_plot(.feature, ncol = 5) |>
  save_plot('sle.skin.dc.inflamm.violin.pdf')

